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INTRODUCTION 


RATIONALE FOR YET ANOTHER ELEMENT 

The discipline known as fracture mechanics has evolved over the last 
half-century from A. A. Griffith’s famous results to the present situation 
of testing standards, advocacy of competing fracture criteria, and potent-- 
and sometimes conflicting- -methods for experimentation and analysis. All 
this work has the common characteristics that first, the presence of a 
sharp macro -crack is presumed, and second, that there is a region of some 
small but finite radius surrounding the crack tip to which we are virtually 
blind. While it is true that various researches have over the years reduced 
the size of this region, the fact remains that we know very little of what 
happens at the crack tip. Of that which we do know, e.g., 'slipping-off’ 
mechanisms, void growth, stress and strain 'singularities’ and so on, it 
has proven quite difficult to translate such information into a more compre- 
hensive and therefore useful form. In particular, our experimental ly -based 
information currently outstrips theoretical insight; this discrepancy prompts 
the present effort. 

What is intended here is to devise means for an improved articulation 
of the analysis of elasto-plastic flow in the vicinity of a crack tip or, 
for that matter, any sharp re-entrant comer. This is to be done through 
the development of a ’special crack-tip element' to be used in conjunction 
with finite element analysis. The element must accomodate those models of 
displacement associated with a crack without at the same time presuming any 
additional features of the behavior which cannot a priori be proven. There 
is thus required a certain economy in the formulation for we are constrained 
to a minimum of input and seek a maximum of output. 



We cannot, for example, be limited solely to elastic behavior as 
characterizes the work of Byskov [1], Wilson [2], among others. Nor may 
we begin with 'fully plastic' behavior as done, say, by Hilton [3] follow- 
ing the work of Hutchinson [4], and Rice and Rosengren [5]. While both 
of these limit cases have their usefulness, the fact that they are limit 
cases precludes information on the process whereby the material goes from 
one to the other. Even the work by Tracey [6], which attempts to include 
both, bypasses the issue of process and blinds us to the rate, means, and 
conditions for the material's transition from purely elastic to fully plas- 
tic response. 

At the present, it appears that sufficient capability and information 
are in hand to proceed along the lines envisaged at the outset of this effort 
A fair amount of experience with and theoretical understanding of the theory 
of elasto-plastic flow has been accumulated so that new means of problem- 
solving may be employed. And, we have enough insight, hopefully, into crack 
behavior to begin practicing the economy outlined above. Some steps in this 
direction have now been taken, and this report outlines the formulation, 
progress to date, and the next steps to be made. 

CONTEXT 

E LAS TO -PLASTICITY IN PLANE STRAIN 

The theory of elasto-plastic flow has been elucidated fully else- 
where (see, e.g., [7]), so that we need only transcribe relations perti- 
nent to our present interest. The theory is stated in terms of increments 
of displacement and stress, and of instantaneous or accumulated values 
of stress. The equilibrium equations , in the absence of body forces, are 
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( 1 ) 


36o /3x + 36t /3y = 0 

x xy 

36x x y/3y + 35o^/3y = 0 

and the constitutive relations for plane strain are written as 
2 u(l + 6s^)«e x = f - v ♦ g(s 2 + 2vs x s z + s 2)]{a x + [-v + e(s x S y - 2vs 2 )]6a y 

+ 2 [ e ( S x + VS z) S xJ 6 V 

2y ( 1 + Ss z) 5e y “[ _v + 6 (Vx ■ 2vs z)] 6a x + l 1 ■ v + K S 5 + + s z JK 

+ 2 [ 6 ( s y + " s z) s xy] 5x xy P) 

4 + 6s z ) 6Y xy = [ es xy( s x + VS z )] 6o x + Ky( s y + VS z )]% 

+ 1 ♦ 2s["s 2 + (1 + v)s 2 l 6 t 

( L xy zj xy 

with the inverse 

6a x /E = £(X + 2>0/E -s 2 /U + 6e x + fx/E - s x s y /(l + 1/6)J 

-[ s xV° + 1/P) ] SY xy 

Sa y /E = ^X/E - Sy s x /Cl ♦ l/6)J«c x + £(X + 2y)/E - s 2 / Cl + l/f!)J<5e y (3) 

-fs s / (1 + 1/0)1 6y 
I y xy J r xy 

5x xy /E = '[ S xy S x /Cl + 1/6) ] Se x '[ S xy S y /(1 + 1/6) ] i£ y 
*[p/ e - s^/d + l/B)j6 Yx y 

In (2) and (3), 6e = 3Su/3x; 6e = 36v/3y; Sy = 3Su/3y + 36v/3x; y is 

x y xy 

the elastic shear modulus; v is Poisson's ratio; E = 2(1 + v)y is Young's 
modulus; X = 2vy/(l - 2v) is Lamp's constant; 0 is the ratio of the elastic 
shear modulus to the slope of the octahedral stress-octahedral plastic strain 
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curve, i.e., 8 = u/u ^ and must be non-negative and bounded; and 
s x = (2o x - o y - a z )/(3^3T o ) 

s y = < 2a y - °Z - ° x )/( 3/3t 0 ) 

s = (2o - a - a )/(3/3 t ) 

z z x y J ' ^ o 

s = x /(/3 t ) 

xy xy v o J 

with 

= (2/9) (a^ +a^+o^-aa - a a - a a +3 t^) 

o x y z yz zx xy xy^ 

Inserting (3) into (1) and replacing fie , 6c , fiy by their equivalents 

x y xy 

in terms of gradients of the displacement increments yields two quasi- 

linear partial differential equations for fiu, Sv. Boundary conditions 

may be written in terms of these displacement increments or the traction 

increments fit , fit which are obtained by using (3) with Cauchy’s relations, 
x y 

Alternatively, one may arrive at a problem statement through use 
of an analogue to the theorem of minimum potential energy. If we define 

n = CE/2)/ [(X + 2y)/E - s 2 /(l + 1/6)] 6c 2 - 2[s y s xy /a + l/6)]s Ey « Y 


[(X + 2 W )/E - s 2 /(l + l/6)]6c y - 2^s x s xy /(l + l/gjJsc^Y 


xy 


xy 


+ [ W/E ' S xy /(1 + 1/g) ] 5Y xy + 2 [ X/E ‘ s x S y /(1 + l/B)]6e x 6 £y |dA 


/ 


(fit 6u + fit fiv) dS 
x y 


where is that portion of the boundary S on which traction increments 
are specified, then an extremum of II (actually a minimum) leads to an 
equivalent result. This approach is germane to any finite element 
formulation, in that the final stiffness equations develop most naturally 
by minimizing the functional H. In the case at hand, this procedure is 
used in a direct manner as described in the next section. 
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DETAILS 


ELEMENT FORMULATION AND PROGRAM STRUCTURE 
At the crack tip, we consider an array of sectors which together 
comprise a special element. Each sector is connected to a regular finite 
element, in the present instance a const ant -strain or linear-displacement 
element. In Figure la, we show a particular grouping of sectors and 
regular elements, and this arrangement has been used in the computations 
described below. In Figure lb, we show details for a typical sector. 

It occupies the space 0 < r < 2r e , 0^ < 9 < 0 2 35 s h° wn ; has nodes at 
( r e > ' (2r g , e i^ * ^ r e J * ( r e * (0 > 6) , usually referred 

to in that order; and has a displacement field given in terms of cartesian 
(u,v) or circular (u^, v ) components. It is convenient always to write 
coordinates in terms of (r, 9) rather than (x, y) . 

The cartesian components of the displacement vector are taken to be 
u = Arcose + Crsine + (E + Fe)r^ cos0 - (G + H6)r^ sine + u 

° (5) 

v = BrsinQ + Drcos 9 + (E + F0)r^ sine + (G + H0)r^ cos0 + v q 

and the corresponding cylindrical components are obtained from 

u = ucos 0 + vsinG 

r 

v = -usin0 + vcos 0 
0 

s o th at 

2.2 

u = Arcos 6 + Brsin 0 + Crsin8cos0 + Drsin0cos0+ 
r 

+ (E + F0)r^ + u cos0 + v sin0 
o o 

2 2 

v„ = -ArsinOcosG + Brsin0cos9 - Crsin 9 + Drcos 0 
u 

+ (G + H0)r^ - u sin0 + v cos0 
o o 

where A, B,..., H are coefficients to be determined; p, q are exponents 
to be found; and u q , v q are rigid translations of the element. In 
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addition, (D - C}/2 is a rigid rotation about the origin in Figure lb. 

In addition, we may note from (5) that A, B give uniform normal strains 
in cartesian coordinates, and that D + C gives a uniform shear. To the 
extent that p, q / integer*, then, terms governed by E, F and G, H rep* 
resent the special aspect of the element. The first pair plus the ex- 
ponent p relate to radial motion, as is evident in (6), and the second 
pair plus q concern circumferential motion. 

Obviously, more complex functions than those shown in (5) or (6) 
may be set up to represent the displacement field within the sector. 

One could add, for example, quadratic and even cubic terms to (5) , or 
the expressions in (5) or (6) may be replaced by isoparametric [8] forms 
following, say, Barsoum [9]. The present objective, however, is not so 
much achievement of high accuracy in the detailed modelling of crack- 
tip behavior; rather our intent is to make sure that a certain approach 
to the problem works at all. Thus, use of a simple displacement field 
seems warranted at the present stage of development. In keeping with 
the economy of formulation noted at the outset, we observe that the ex- 
pressions in (5) or (6) possess the minimal attributes of rigid motion, 
uniform strain, and the special feature we seek to establish [8]. 

The next step in formulation is to replace the coefficients A, B, 

. . . , H in (5) and (6) by their equivalent in terms of nodal point dis- 
placements. We therefore construct a vector of nodal displacements. 


referred to the cartesian coordinate directions, as 
{u T } = {u x v : u 2 v 2 u 3 v 3 u 4 v 4 u 5 v s > 


★ 

Note that 0 < p,q < 1 is the operating range of interest. 
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In (7), superscript T denotes transpose. The order of nodal numbering 

begins, as indicated above, at (r , 8,) and continues counterclockwise 

e j. 

to r = 0. Then, using the vector 

{A 1 } = {ABCDEFGHu v} (8) 

— o o 

we may write (5) at each of the five nodal positions to arrive at a set 
of ten equations of the form 

{u} = [Q] {A} (9) 


where [Q] is a 10 * 10 matrix whose elements depend on nodal point coordi- 
nates. Inversion of (9) leads to 


{A} = [Z] {u} 

Then (5) may be replaced by 

“ £8)- <h} 


and (6) by 


u r (r,9) 

v Q (r,e) 


= [Y (T j 0) ] (u) 


( 10 ) 


(5’)- 


( 6 ? ) 


where [a(r,0)j and [y(r,0)] are 2 x 10 matrices whose elements 
nodal coordinates and position of reference within the sector. 
we may differentiate (S’) to obtain a vector of strains, viz. 

= [B(r,e)3 (u) 



depend upon 
Finally, 


(ID 


* 

and [3(r,8)] is a 3 x 10 matrix. Elements in the pertinent matrices 
listed above appear in the Appendix. 

In the same manner that the strain- displacement relations are trans- 
formed from total or accumulated values, to incremental values, we may 


The matrix [8] in (11) is distinct from the modulus ratio in (2) and (3). 
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alter (S’) and (11) for use in elasto-plasticity . Thus, if the incre- 
ments of displacement at nodal points may be formed into a vector (dp.) 
in the pattern of (7), the displacement increments interior to a sector 


are given by 


6u (r,8) 
6v (r,0) 


[a(r , 0) ] {<$u} 


Furthermore, the strain increments are written as 

- [0(r,8)] {6u} 


'6c x (r,0) } 
{6e} =^c^(r,0) 

,'SY X y.(t,0) 


( 12 ) 


(13) 


and the matrices [a(r,0)] and [0(r,6)] are the same as given in (5 1 ) and 
(11), and as detailed in the Appendix. 


We introduce next a vector of stress increments as 
[6o x (r,9) 


{6a} ={5a^(r,0) 

^T xy (r,6), 


(14) 


and, from (3), it follows that • 

{6a} = [M] {6e } (15) 

where [M] is a symmetric 3x3 matrix whose elements depend upon the 
accumulated or instantaneous stresses, the elastic constants, and slope 
of the octahedral stress -octahedral plastic strain curve. Association 
of (3) and (15) indicates that 

Mjj = A + 2v - E/Cl + 1/8) 

M 22 = A + 2p - s 2 E/Cl + 1/8) 

M 33 >* - 4 E/ a + 


“ ^21 ^ s x S y ^/ (1 + 1 / 8 ) 

M 23 = M 32 = ■ Vxy E/C1 + 1/S) 

M 31 " M 13 = - Vxy E/(1 + W6) 
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and the definitions of the various terms remain as stated in the pre- 
vious section. 

The functional in the integral theorem cited above now takes the form, 
for one sector. 


®2 2r e 


= .j f J {Seb [M] {Se} rdrdS - {if} {{u} 


Q x 0 


or 


0 2 2r e 


= J {fiu} // [e T (r,e)] [M] [g(r,9)] rdrdB {6u} - {«f 


e o 


{6T‘> (6u) (16) 


in which {<57^} is a vector of nodal forces in an order that corresponds 
to that used in constructing {5_u}--see (7). 

At this stage, the normal next step is to require that IT be station- 
ary , with respect to nodal displacement increments, and the result is a 
linear relation between {6u} and {5T} viz, 

[K g ] {6u} ~ {6T} ' (17) 

where, for the sector, 

2r 

2 e 

I K S ] •// [3 T (r,6)] [M] [B (r , 0) ] rdrde 


While. the same sort of result obtains in the present case, two additional 
relations must be developed. That is, we must determine conditions in 
p and q, as used in (5) or (6), from suitable manipulation of n. Let us 
suppose, for example, that each sector has its own values of p and q. 
Then, in addition to (17), we must require 

3H /3p =0, 3IT /3q = 0 
s s 
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It may be shown, however, that this supposition leads to incompatibilities 
between adjacent sectors of the special element, to the extent that p and 
q vary from one sector to the next. Hence, in a formal sense, it is pre- 
ferable to define 
n 

"e -J/s U8) 


where n is the total number of sectors in the special element. Then, 
letting p and q be common to all sectors, we require that 

5IE e /3p = 0, 3n e /3q = 0 (19) 

which gives the additional relations needed. 


The foregoing gives a useful clue to tactics for problem solving, 
and a variety of alternatives may be listed. Let us suppose that ex- 
pressions of the form (16) are summed via (18); to that total, we add 
values of for all m regular elements surrounding the special element, 

viz: 

n ra ‘ ' 


n = £ n - 

S= 1 S 



r=l 


( 20 ) 


The final value H is evidently quadratic in the nodal displacement incre- 
ments throughout the domain modelled for a given problem and horribly non- 
quadratic in the exponents p and q in the special element. Because the 
overall problem is thus non-quadratic, one could approach its solution 
by attempting to minimize JI (in (20)) directly. This has the advantage' 
of being notionally straightforward; it is at the same time quite expen- 
sive because algorithms for achieving extremal values are operationally 
slow when large numbers of variables are being treated. 

Alternatively, we may utilize the fact that most of the variables 
in (20) appear in quadratic form. Writing the symbolic expression 
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dH = (3JI/3{6u}) d{<5u) + (3II/3p)dp •+ (3ll/3q)dq = 0 
where {<$u} (in contrast to {611}) is the set of all nodal displacement in- 
crements in the domain of interest, we observe that a sufficient condition 
for the solution is 

3H/3{<Su} = [K] (6u) - {6T} = 0 

an/3p = sn /3p = o ( 21 ) 

3ir/3q = an /aq = o 

© 

In (21), [K] is the master stiffness matrix for the domain, {6T} is the 
set of all nodal force increments, and the first relation is thereby a 
set of linear algebraic relations for (Su) in terms of {6T}, the latter 
presumably known. The second and third relations, however, are each single 
non-linear algebraic equations governing p and q. All three of (21) 
must be satisfied simultaneously. 

Using this alternative approach, we have .constructed a computer pro-" 
gram SPECEL for the special element embedded in an array of regular el- 
ements. Its basic procedure for the kth increment may be outlined simply 
as follows. The increment begins by setting p and q to the values ob- 
tained from the previous increment and, using these values, a set of stiff- 
ness equations is set up. Thus the first of (21) is solved. The energy- 
like term is found and minimized with respect to p and q, that is, the 
exponents are adjusted to effect the second and third of (21). With these 
new values of p and q the stiffness equations are re-solved. This pro- 
cess is repeated until a minimum of n itself is established so that (21) 
is satisfied. The results are then processed in the usual manner [10], 
checking to see that any unloading is taken properly into account, com- 
puting and printing various items of information such as displacements, 
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stresses, and the J-integral. A flow chart is given in Figure 2. 

At the present writing, an exception to the foregoing procedure is 

programmed for the first load increment. Since this step is wholly elastic, 

the exponents are pre-set to precisely one-half (0.5) and no searching is 

performed to find 'better 1 values, i.e., values which would reduce II . 

e 

This exception stems from efforts to verify performance of the code by 
comparing it to known elastic solutions; we have yet to decide whether to 
retain this feature or, through a trivial change in code, force the first 
increment into the same pattern as obtains for subsequent load steps. 

STATUS 

ELASTIC AND INITIAL PLASTIC RESPONSE 

Having written code to implement the formulation described above, 
the next step was to check it out. in'isome detail. It' happens that, 
although the bulk of the formulation involves the concepts expressed in 
the preceeding section, the bulk of the code itself involves the details 
of matrices such as those given in the Appendix and their manipulation. 

Thus verification can and did go forward in terms of elastic analysis, 
i.e., using the constant terms in [M] given in (15) et seq f 

Using both the element arrays shown in Figures la and 3, a sequence 
of problems was solved. The first set involved uniform tension and 
uniform extension of the two maps. Results of these analyses indicated 
some errors in the code which were then corrected. What was determined 
ultimately is that the special feature of the element creates a low noise 
level in the sector's force vector; the singular terms are active, in 
effect, even though excitation is solely uniaxial homogeneous stress. 
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Figure 2 - Flow chart for SPECEL 

Read input data; initialize parameters, arrays 
INCR = 0 

10 INCR = INCR + 1 
I CYC = 0 

20 I CYC = I CYC + 1 
ITER - 0 

30 ITER = ITER + 1 

Generate stiffness matrix [K] 

Apply boundary conditions {ST} 

Solve equation set {Su} 

Obtain strain increments {6e} = [8]{<5u} 



r~- 




Figure 3 - Special element embedded in a larger array of regular elements 



I 


The magnitude of resultant noise is, however, quite low. If, for example, 

the applied stress is a and or = 28125, the noise level nowhere exceeds 0.32% 

e 

of this amount when 24 sectors are used in the special element. This effect 
was examined in some detail. We find that it is predictable in terms of 
the formulation (5), that it is of the order of A 2 (where a = 6 2 - 6^, and 
that it is not an objectionable result of the simple representation (5) 
of a complex situation. 

Moving next to crack problems, the map in Figure 3 was then subjected 
to two loadings. The first is uniaxial tension of a cracked panel, where 
the crack length/ ligament length ratio is 1/2. For that case, the con- 
ventional stress intensity is given [11] as 
K = 3.17 a/a 

where a is the crack length, here set to unity. The computation gave an 
average value of J over seven paths, which when converted to Kj , is 
Kj = 3.122 d/a 

and the coefficient of variation of Kj is 0.0090. Thus the error in stress 
intensity is about 1.5% and the variation* is about 0.9% over circular 
paths having radii ranging from 0.0032 to 0.1218. (It may be noted that 
crack length is 1 and that the special element radius is 0.125.) This 
result is regarded as satisfactory. 

The second loading imposed on the range of Figure 3 may be termed 
a ’pure’ Kj excitation. At all boundary nodes except those along the 
flank of the crack, it was required that 


The sample coefficient of variation is the standard deviation divided by 
the average'. 
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( 22 ) 


u = /(2r) x 10 ^ x (5/3 - cos 0 ) cos 0/2 
v = /( 2 r) x 10 2 x ( 5/3 - cos 6 ) sin 0/2 

The values above result if E = 10^ lb/in^, v = 1/3, K^. = I.s/ft x 10^ lb/in^ 2 , 
and the problem is in plane strain. Using the procedure outlined above 
to obtain K^., a value of 1.5363/ft * 10 Ib/in ■ was achieved, the coefficient 
of variation being 0.92%* The error in is seen to be 2.42% which, we 
believe, is explicable as follows. A complex (i.e., non-homogeneous) ex- 
citation is transmitted from the outer boundary to the special element via 
regular elements. The latter are known to be slightly stiffer than the 
continuum they represent so that, in the present instance, the loading * 
on the special element will exceed the precise amount. We observe here 
that the excess is modest, and we expect that this pattern will be reflec- 
ted in the computed nodal point displacements relative to the exact solu- 
tion as given in (22). As seen in Figures 4 and 5, such is the case. . 
Moreover, we observe that the circumferential behavior of the displacements 
is quite accurate and, in terms of nodal displacements for the special 
element, response is quite satis facto ly. 

Stresses interior to the element are also acceptable. In Figures 
6 and 7 we plot the principal stresses which may be derived from (22) as 
a.. = 1.5 x 10^ x /(a/2r) x (1 + sin0/2) cos9/2 

(23) 

- 1.5 x 10 ^ x /(a/ 2 r) x (1 - sin 8 / 2 ) cos 9/2 

Values derived from (23) are given in Figures 6 and 7 for three values of 
r/a, along the midline of each sector. It is seen that the stresses tend 
to be high but not by an amount that greatly exceeds the original over- 
loading through amplification of . Moreover, the circumferential grad- 
ient is well replicated except where the stresses approach null values 


- 17 - 



/ / (2 ar) x 10 


Figure 4 - Horizontal displacement, normalized, for 

pure Kj loading at remote boundary as a 

function of angular position- -nodal point 

values. Kj = 1.5/ir x 10 5 lb/in^ 2 , 

r / a = 0.06250, v = 1 / 3 , E = 10 7 lb/in 2 , 

© 

plane strain. 

a nodes at r = 2r (2 7-51) 
e 

o nodes at r = ( 1-25) 
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01 X {xi ?z)f/ 


2.5 h 


2.0 r 


Figure 5 - Vertical displacement, 

normalized, for pure loading 

at remote boundary as a function 

of angular position--nodal point 

values. Kj = 1.5 /tt x 10 5 lb/in 3 / 2 

r /a = 0.06250, v = 1/3, 
e 7.2 

E = 10 lb/in , plane strain. 


O £ 

D * # 


1.5 h 



— j- — 


1.0 


a nodes at r = 2r & (27-51) 

B nodes at r = r ( 1-25) 

6 

• exact solution v^//(2ar) = 

10' 2 (5/3-cos 9) sine/2 


0.5 


30 


60 


90 


120 


150 


180 


9, deg 
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/(2r/a3/ Y 


Figure 6 - Larger in-plane principal stress , normalized, 
for pure K T loading at remote boundary as a function of 
angular position-element values. = 1.5 /jt x 10 lb/in 
r /a = 0.06250, Y = 10 5 lb/in 2 . 





A r/a = 0.087865 
□ r/a - 0.037135 


v r/a = 0.003181 


• exact solution a^/(2r/a)/y = 
1.5(l+sine/2) cos 0/2 


8, deg 


A/ (*/*Z)/ 


Figure 7 - Smaller in-plane principal stress, normalized, for 
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along the flank of the crack. 

With indices of performance such as those shown, it is concluded that 
the elastic response of the program SPECEL is satisfactory. Improvements 
are accessible through reduction in sector angle, i.e., increase in number 
of sectors, or through refinements in the interpolation functions (5) or 
(6). Such i up rove men t , while eventually worth seeking, is separate from 
the matter of central interest at this stage, namely, how well the overall 
scheme works. 

As an initial foray into this next stage, we adjust the initial 
loading so that yield is just beyond the point of initiation at the most 
highly stressed point in the body. For uniform tensile loading this 
occurs in the twelfth sector and the radial location of the point is 
r/ a = 0.003181. The next load increment is set at 5% of the first, the 
third at 5% of the sum of the first two, and so on. In this manner, 
loading is built up of a succession of increments [7,10], and we have 
looked at results for a total of five loading increments. 

We observed that both exponents decreased slightly as loading 
progressed, from the pre-set initial values of 0.5 to p = 0.402 and 
q = 0.442. It is reasonable to expect that the exponents should decrease 
since that implies an increase in the singular nature of the strain in- 
crements which, after all, is the sort of behavior accompanying yield. 

It should be borne in mind that these values of the exponent refer to the 
displacement increments associated with the fifth load step, not to an 
accumulated value of the respective exponents. We anticipate that, as 
yield begins, the exponents will change rapidly as the material accomodates 
the changeover from purely elastic to elasto-plastic response; thereafter 
the rate at which the exponents change should decline. This pattern of 
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response characterizes other plane strain analyses we have observed and is 
discussed at some length in [12], 

The second observation made from the initial foray into elasto-plastic 
response is that the computational time required per load increment is 
greater than we prefer, about 1 min on the CDC 6600 and over 2 min on the 
Univac 1108. While such, amounts of time are not prohibitive, they are 
uneconomical and steps are now being taken to reduce the computer require- 
ments for the program to operate. 

CONCLUDING REMARKS; A CAVEAT 

The formulation of a ’special element 1 for analysis of elasto-plastic 
behavior in the vicinity of a crack tip has been outlined, and progress 
to date summarized. Basically it is evident that the approach is oper- 
ationally feasible, although certain improvements to code are required 
before the present program is attractive. Once done, however, the pro- 
cedure should permit an examination of the process of yield in the crack- 
tip vicinity, taking properly into account the continued presence of 
elastic strains, plastic strains according to an arbitrary work -hardening 
stress-strain relationship, unloading, and the sequential nature of the 
flow process. 

It has been noted [13] that an aspect of this approach is open to 
question, in that the information sought should overlap that already given 
by the HRR model [4,5], as it is usually termed. The commentary appears 
to take the following form. If the element size is large with respect to 
the plastic zone, then the result would be essentially an elastic singularity. 
Conversely, if the plastic zone is large with respect to the element, the 
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HRR result should obtain in that it is an asymptotic solution, i.e., it 
describes behavior within small radial distances from the crack tip. An 
intermediate situation would perforce give special element results that 
are some sort of average of the two types of solution. Viewed in another 
way, this comment suggests that, since the special element is of fixed 
size in any given analysis, its result in terms, say, of p and q should at 
most go through a smooth transition from purely elastic values (1/2) to 
those associated with the HRR results. Hence, there is little to be 
gained by pursuing this formulation. 

Some response to this commentary is in order. There are two aspects 
to the problem under consideration. One is the (changing) structure of 
the singularity, and the other is the rate of change it undergoes. The 
structure is initially that of elasticity and is long since established 
[14]. It then is altered to something else which, for one somewhat special 
set of conditions, is also documented [4,5]. Not known, however, is how 
that second state may be affected by the presence of elastic strains in 
addition to plastic strains, alternate material behavior (in the sense of 
the stress-strain curve), and use of flow theory in place of deformation 
theory as a model of the material. It may be that these factors have but 
a marginal influence, or they may affect the result to a more pronounced 
degree. Putting the matter in such simplistic terms, however, begs the 
second point which may be viewed as concern with the structure of the sin- 
gularity at any of a sequence of 'intermediate’ loads. It is not clear, 
for example, that the singularity goes monotonically from purely elastic 
response to 'fully plastic’, that there is not some alternation associated 
with two or three phases of load redistribution as yield and therefore 
' hardening progress. To the extent that such response occurs, it is natural 
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to expect that both intensity and structure of a singularity will be sensi- 
tive to it; none of the information we have seen in the literature reflects 
this sort of sequencing of events, or process. In spite of the commentary 
[13], the present work appears worth pursuit, although we stand reminded 
that sizing the special element is a sensitive matter which bears close 
examination once we reach the stage of making analyses on a productive basis. 
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APPENDIX 


ELEMENTS IN MATRICES 


Using the abbreviations 
= sin0^ 

s 2 = sine 2 


= COS0^ 

c 2 = cose 2 


The non-zero elements of Q are 


*11 

^24 = 

r c 1 
e I 

Q 22 

" Q 13 = 

r s 
e I 

*31 

Q 44 

2r C- 
e 1 

Q 42 

= Q 33 = 

2r e s l 

*S1 

= Q 64 = ' 

2r c„ 
e 2 

%2 

= ^53 = : 

2r e S 2 

Q 71 

= Q 84 = 

r e C 2 

^82 

= Q 73 = 

r s 
e 2 

*is 

= 

- r^c 

e 1 

q 17 

= Vi 

q 

= -rs 

e 

^25 

= W 0 ! 

= A, 

. e I 

Q 2 7 

" <W 9 1 

= r^c- 

e 1 

^35 

= <V e i 

- 2 p r P c 
e 1 

Q 37 

= Q 38^ 6 1 

q 

= -2 n r s. 
e 1 

Q45 

Q 46^ 6 1 

= 2*1*3, 
e 1 

Q 47 

= W 0 ! 

0 q q 

= 2 r e c l 

Q 55 

= <V 9 2 . 

- 2 ^2 

%7 

= Q 58 /6 2 

= -2 q r q s 
e 2 

^65 

= <W e 2 

= 2Pr ^ S 2 

%7 

= <W 0 2 

= 2^r^c 

e c 2 

Q75 

= Q 76 /0 2 

" T e c 2 

Q 77 

= Q 78 /6 2 

q 

— -r s 

e 2 

Q 85 

= <V 0 2 

= x*s 

e 2 

Q 87 

= ^88^2 

q 

= r n c^ 

e 2 

V 

9 = Q i+1, 

= 1 

10 * 

i = 1 

3 5 7 9 
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More abbreviations are useful in detailing the elements of [Z] 
they are 


P 2 = 1/(2 - 2 P ) 
Q 2 = 1/(2 - 2 q ) 
P = P 2 /(r^A) 

Q Q 2 /(A) 


R = 1/Cr^sinA) 

R r " C P 2 - V R 

R 1 = R r s l s 2 c l 
**2 = R r S l 5 2 C 2 


R = R see 
3 r 1 1 2 


R 4 = R r S 2 °l C 2 


G k = t p 2 c k + Vk> R 

H k - (p 2 s k + 



1,2 


The non- zero elements of [Z] are given by 


Z,, - Rs~ 
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1 2 



- 2 Z. _ 



2 R_ 

21 


23 


3 
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2 G, 

31 2 
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1 2 
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- 2 Z . «. 



- 2 R, 

41 


43 


1 


_ 

- 2 Z_ _ 

_ 

2 Pe_c. 

51 


53 


r-t 

z„. 


- 2 Z. , 



- 2 Pc, 

61 


63 


1 

Z 71 

= 

- 2Z 73 

= 

- 2 Q 8 2Sl 

Z 81 

= 

- 2Z 83 

= 

2 Qs j 


-2Z, r 
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- 

2G_s _ 

15 


17 1 


2 1 

-2Z„_ 




-2R . 

25 


27 


4 

-2Z_ r 


Z — Re, 

_ 

-2G_ c. 

35 


37 1 


r-i 

rN. 

i 

i 

-2Z ir 



Z , „ 

_ 

2R„ 

45 


47 


2 

-2Z__ 

_ 

Z rn 


-2P8, c„ 

55 


57 


1 2 

- 2Z 6S 

= 

Z 67 

= 

2Pc 2 

- 2Z 75 

= 

Z 77 


2Q0^S2 


' 2Z 85 = Z 87 = ‘ 2Qs 2 
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Then 


12 


= - 2 Z = - 2 R- 

14 1 


2Z 16 Z 18 


= 2 R. 


Z 22 + RC 2 = 


-21 


24 


2 H 1 C 2 


- 2 Z 


26 


Z 28 " Rc l 


“ 2 H 2 C 1 


'32 


' *2Z 34 = 2R : 


“ 2Z 36 = Z 3S 


= - 2 R 
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- 2 Z 
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'52 


= -21 


54 


2 P 6 2 S 1 
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- 2P9 1 S 2 


'62 


= - 2 Z . . = - 2 Ps. 
64 


- 2Z 66 = Z 68 


= 2 Ps, 


72 


'82 


= -2Z 74 = 2Q0 2Cl 
= - 2 Z g4 = - 2 Q Cl 


' 2Z 76 “ Z 78 
’ 2Z 86 = Z 88 


■ -2Q0 x c 2 
2 Qc 0 


Z j ,9 Z j, 2 n -1 

n=l 

Z j , 10 “ ^ - Z m, 2 n 

n - 1 


j = 1 , 2,. . . ,8 


7 =7 =1 

9,9 10,10 


The interpolation functions [a(r, 0)3 may be written in terms of 


P r = PrP = (P 2 /A) (r/r e ) P 
Q r = Q T<1 = (Q 2 /A) (r/rj q 


d 2 = e 2 . - e 


d i ■ 9 - 9 i 


s = sin 0 
c = cos 0 


S 2 = Sind 2 


« sind^ 
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21 
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a 12 

= 

- 2a 14 = 
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1,9 = 1 
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4 
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Next, let 

fjCe) 


V e > 


d^(s 2 + pc 2 ) - sc 
l 2 


2 2 
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P » 
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P /r 
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So that elements of [0(r,8)] become 
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Finally, we note 
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is the simple transformation required. 
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